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Abstract — Heat-transport is important for geothermal exploration. The presence of 
fractures can have a pronounced effect on groundwater and heat transfer .The inclusion of 
fractures into geothermal reservoir models on different scales is often still a difficult task. A 
comparison of approaches for flow in fractures has been carried out. A very simple 
approach is to simulate fractures with thin but highly conductive layers, for instance by 
applying the Cubic-Law. A more sophisticated approach, typically in FEM codes, is the 
application of lower dimensional (1D/2D) high permeable discrete elements with specific 
flow properties, following e.g. Hagen-Poiseuille or Manning- Strickler. However, such an 
approach typically fails while studying only partly saturated fractures. For studying the 
applicability of simplified fracture modellingapproaches a comparison with a CFD 
(Computational Fluid Dynamics) solution was performed. Furthermore a DEM (Discrete 
Element Method)approach has been illuminated. The various methodologies are studied by 
varying roughnesses,this way studying the versatility of the approach. The sensitivity of 
flow in fractures to various numerical parameters can be studied this way. A detailed 
analysis of temperature and flow using Peclet and Reynolds numbers helps to quantify the 
contributions of the different transfer processes. 

Index Terms — Fractures, Cubic Law, Discrete Elements, Computational Fluid Dynamics, 
Finite Element Method, Discrete Element Method 



I. Introduction 

Geothermal energy is a cost effective, renewable way to process transfer energy and an important technology 
with respect to a change towards renewable energy. The increasing power demand, especially in developing 
countries,together with the limited availability of fossil fuels has led to increased effort in this field. The 
basic idea geothermal energy is to exploit the heat generated in the depth of the earth's crust, mainly by 
natural radioactive decay. The heat can be extracted by the circulation water in greater depths. Outside of 
volcanic regions a depth of more than 3 kilometers is necessary. It involves the drilling of boreholes, 
followed by the pumping of water to depths where the temperature of the rock is markedly higher than on the 
surface. The resulting vapor and water is collected by another borehole and passed through a steam turbine to 
produce electrical energy. The presence of fractures has a strong impact on the subsurface temperature 
distribution [1], 

A more recent analysis compares the advances of geothermal modelling providing a holistic view on the 
current state of research involved [2]. Modelling of geothermal reservoirs is possible with a number of 
different codes. In this study, the use of the software FEFLOW 6.1, which is software devised by DHI 
WASY,has been proposed which will be used to analyse the situations of heat transfer and groundwater 
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flow. FEFLOW allows the modelling of heat transport in the subsurface and the inclusion of fractures using 
discrete feature elements. The software and its capabilities have been illustrated multiple times by Ref. [3] 
and [4] . The modelling to simulate the actual network of fractures would be ideally done using the cubic law 
which will link the velocity to the fracture aperture and pressure. Other methods were employed due to the 
lack of computational power. 

The modelling of flow through fractures bears an important role in the geothermal industry. With the advent 
of computational power the reservoir simulations in the geothermal industry are very useful in both in terms 
of industrial production as well as scientific study. The FEM codes like FEFLOW help in the simulation of 
heat transport and flow in a given geothermal reservoir. However in these simulations it is important to 
simulate the flow through fractures systems as fractures reservoirs encapsulate an imperative part of 
production and their geological study has been explain in detail by [5]. There is a marked difference in the 
methodology of production and extraction from a fractured reservoir which has been dealt by Ref. [6] . Hence 
from a scientific perspective it remains important to look into the modelling of flow in reservoirs as the flow 
in a reservoir rock cardinally controls the production. A stochastic numerical model was built and the 
flowtransport in fractured rock was penned out by Bear Ref. [7]. There have been numerous methods 
developed over the years to model the flow of fluid and heat transport in fractured rock by Ref. [8]. However 
in FEM codes like FEFLOW, the introduction of fractures is through discrete elements in the model of the 
reservoir. The modelling of the flow in these discrete elements can take various paths as the flow law can be 
varies according to the nature of the fracture involved. The finite element code has been elucidated and its 
applications for fluid flow in rock mases described by Ref. [9]. The issues and examples of geothermal 
processes has been described by Ref. [10]. However this maybe computationally intensive and the use of 
discrete elements is not possible in the case if unsaturated flow simulations. Hence an alternative is to look at 
the applications of CFD which may be used to model the flow in fractures. This can be used to model the 
flow in the fractures, serving to be both computationally efficient as well as accurate. 

Analyzing results can be a complicated task, which has to be maneuvered in the direction of the field of 
applicability. Hence in this paper we will look at the temperature plots and the hydraulic head plots of the 
simulations to predict the outcome. Another interesting observation can be drawn using the Peclet number as 
shown by Ref. [11]. Since the validation of the flow in fractures cannot be carried out without the inclusion 
of field data, an approach through computational fluid dynamics has been sought. Computational Fluid 
Dynamics (CFD) is the art of simulating a flow of fluid in a given geometry. It is often used to assess the 
performance of engineering devices, to explore in a cost-effective manner several competing designs, or to 
provide understanding of flow processes within or around a given configuration. 

The OpenFOAM® (Open Field Operation and Manipulation) Ref. [16] is a CFD Toolbox is afree, open 
sourceCFD software package produced byOpenCFD Ltd. It has an extensive range of features to solve 
anything from complex fluid flows and can be efficiently used in this case to model the flow in fractures. The 
application of CFD to the modelling of flow in fluids has been illuminated in the recent past [12]. The 
modelling of fractures in rock masses has been taken on extensively using Discrete Element Method. It is 
possible to model their yield strength, the shear strength of fractures as well as the cross sectional area of 
fractures in rock using this method. Hence this paper aims to look at the prospects of the software 3DEC - an 
advanced code for the numerical modelling geotechnical analysis in discrete element blocks - to model the 
flow of fluids in fractures for elucidation of modelling flow with the aforementioned simulations in 
OpenFOAM and FEFLOW. 

II. Model Setup 

A. Conceptual Model 

As a first step in model design FEFLOW requires the initial mesh extent where we define the range of the 
model in terms of distance in which the model exists. The model locale or dimensions have been defined in a 
coordinate system of 1 m times 1 m.The following simulations were composed of standard equation for 
saturated flow, and shock capturing was used in the upwinding. 

B. Numerical Model 

The groundwater flow equation is the basis upon which the flow phenomenon is simulated. It is common 
practice to develop numerical solutions of the equation based on the specific boundary conditions. In the 
present case, a finite element grid was defined. The grid provides for the discretization of the data - which is 
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performed in both the horizontal and vertical directions as shown in figure 1. In the present case, the 
horizontal discretization was specified by defining a triangular mesh generator [13]. 




Figurel. Mesh of the Model (Top) three-dimensional model view (Bottom) 



It is necessary to define a set of boundary conditions (BC) for the simulation of the model. By default, all 
model boundaries in FEFLOW are impervious. FEFLOW supports four basic types of boundary conditions 
for flow, mass and heat transport. For the current model we adopt the Dirichlet type BC which requires us to 
define a boundary head condition as well as a boundary condition for the temperature. The temperature BC is 
defined at the injection well as well giving it a temperature of one degree Celsius which acts as a tracer. 
Using groundwater a tracer has been discussed in detail before Ref. [14]. The boundary condition for the well 
as also defined where the pumping rate of the injection and extraction well maybe assigned. The software 
FEFLOW allows the user to define initial conditions for the model in question. Hence the initial conditions 
such as the temperature, hydraulic head, conductivity, porosity, etc. 

III. Simulations 

The simulation is carried out by FEFLOW according to the defined time period after the entire initial and 
boundary conditions are fulfilled. The corresponding charts of temperature, head, pressure, etc. are plotted 
automatically by FEFLOW.The study has been divided into distinct cases based on the flow in fractures 
using discrete features and a highly conductive layer. 

The models are mainly of four types and their flow has been studied below. The well has been given a 
boundary condition with a temperature which acts as a tracer. The models have a pumping and extraction 
well which have been placed 0.5 m apart. The wells are being pumped at a rate of -5 m 3 (injection) and 5 
m 3 (extraction) per day. Fractures in the form of highly conductive slices /discrete features were added in the 
model in the middle slice. The model adopted the Dirichlet type BC which requires us to define a 
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boundaryhead condition as well as a boundary condition for the temperature. The boundary conditions for 
head and temperature was uniformly set to zero across the model. The figures below represent the flow 
results for the models described above. Figure 2 is the model where a highly conductive layer of hydraulic 
conductivity of 100 m/s has been used in the slice for the simulation of fracture flow. Figure 3 represents 
fracture flow using discrete features and apply Manning Strickler law to the flow in the fractures, where the 
aperture used is 0.005 m. Figure 4 represents fracture flow using discrete features and apply Hagen-Poiseuille 
flow to the fractures, where the aperture used is 0.005 m. Figure 5 represents fracture flow using discrete 
features and applies Darcy law, where the aperture used is 0.005 m. The flow equations used in the 
simulations are: 
Darcy flow law: 

q = -~(Vp) (1) 

k is the hydraulic conductivity, (j, is the viscosity 



Hagen-Poiseuille flow law(Cubic form) : 



b is the aperture, (a. is the viscosity. 



12/i dx 



(2) 



Manning-Strickler flow law : 



VllVM ; 



(3) 



h is hydraulic head, a is a constant, r is the aperture, I is the flow gradient and t is the friction factor. 
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Figure 2. Hydraulic head and temperature results for model with highly conductive layer as fracture 

The Peclet number was calculated using an IFM Plugin to the software. The plugin was coded using 
Microsoft Visual Studio and programmed in C++. Pecletnumber was calculated according to the formula: 

Pe (5) 
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A 

a = — 

pc 



(6) 



The characteristic length was computed by calculating the value of the volume of the element divided by the 
area of the top face of the element. This provides an equivalent of the characteristic length in the equation 
mentioned above. 
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IV. Openfoam Modelling 

Groundwater flow simulations based on Darcy's law are typically limited to laminar flow. However in cases 
where geothermal reservoirs are employed the flow can become turbulent in cases of involvement of hot dry 
rock, fractures and geothermal wells where water is pumped at a high velocity. Hence it becomes vital to 
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Figure 5. Hydraulic head and temperature results for model with discrete features using Darcy flow 

examine the flow both in the cases of laminar as well as turbulent flow in OpenFoam. Thus in the following 
models have been created to study the flow using laminar as well as turbulent models. The models have been 
meshed and computed in OpenFoam. The post processing has been performed using ParaFoamtogether with 
ParaView. 

Since OpenFoam is a 3D solver the meshing has been done for a 3D geometry however the analysis itself is 
in 2D as two of the faces have been declared empty. The geometry of the model consists of a rectangle with a 
length to diameter ratio of 15. The length of the model is much larger as compared to the diameter so as to 
make sure that the flow is fully developed in the denoted simulation time. The mesh grading can be specified 
in OpenFoam. The mesh grading was simple and uniform throughout the model. Refinement of the mesh was 
not necessary near the walls as OpenFoam provides the wall function which is an alternative to the 
refinement of the mesh. The meshing of the model as well as a screen of the model has been displayed below 
in figure 6. 



Mesh Information 

boundingBox: (0 8 a) (1.5 9.1 0.01) 
nPoints: 49282 
nCells: 24000 
nFaces: 96640 
nlnternalFaces : 47360 

Patches 

patch (start 

patch 1 (start 

patch 2 (start 

patch 3 (start 

patch 4 (start 



Figure 6. Meshing information 



: 47360 size: 40) name: Inlet 

: 47400 size: 40) name: outlet 

: 47440 size: 600) name: upperWall 

: 48040 size: 600) name: lowerWall 

: 48640 size: 48000) name: f rontAndEack 
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OpenFoam requires the selection of solvers at the early preprocessing stage. The solver used in all these 
simulations is called pisoFoam. PisoFoam has been defined in the OpenFoam libraries as standard solver 
which is used for incompressible flow. It can be used to simulate both laminar as well as turbulent flow. The 
turbulent flow has been simulated using the RAS (Reynolds -Averaged Simulation) package.. 
The models have been built on the pisoFoam which consists of an inlet, outlet, top wall and bottom wall. The 
inlet velocity has been defined, and the inlet pressure in variable. The outlet pressure of the models is zero. 
This thereby permits the user to measure the pressure difference at points of interest which when divided by 
the velocity will give us the required result. 

The model has been run for a period of 10 seconds. Figure 7 shows the plot of the pressure at the end time of 
the simulation. Figure 8 gives the plot of the velocity at the end time of the simulation. These simulations can 
be used for comparison with the FEFLOW simulations of fracture flow. The software uses the K-epsilon 
model for the simulation with the general turbulence momentum and energy equations which have been 
described in detail by Ref. [15]. 
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Figure 7. Pressure at end time with Surface Roughness (5) & Turbulent Flow 
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Figure 8. Velocity at end time with Surface Roughness (5) & Turbulent Flow 



V. Feflow Comparison With Openfoam 

A simple FEFLOW model was built in which flow was simulated without the use of wells or temperature 
conditions. The model was simulated with a flow which is driven by a pressure gradient equivalent to the 
result from the CFD simulation in OpenFoam. The simple groundwater equation has been employed without 
mass or heat coupling. The pressure gradient was fixed to 0.31 kPa and the flow was simulated in four 
regimes namely flow in a highly conductive layer which behaves close to flow in fracture, discrete features 
used to describe fractures using Hagen-Poiseuille, Manning-Strickler and Darcy flow laws. The ratio of the 
velocity and the pressure gradient has been compared to the OpenFoam results. Table 1 and Figure 9 display 
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the nodal Darcy velocities for the models generated in OpenFoam. In the case of the model where a highly 
conductive layer has been used to depict the fractures, it can be seen that there is a uniform flow in the 
middle and center of the model however at the edges there seem to be an increase in the Darcy velocity. This 
is because the boundary conditions have been applied to the borders of the model which results in certain 
edge effects. However this does not influence out study as we are interested in the velocity in the center of 
the model which is not affected by developing flow or boundary conditions. 

Darcy flu* (nodal) 
- Fnnpe* - 
pttj 

■ 30 0928 31 2TS7 

■ 11 mi „ 30 0326 
' 25 m V 7778 



23.1491 
20.3333 
I 18.5165 
l-j ItiV 



15.493 
33.1 4S1 
20 6333 

ia.5ies 



I 113I3B9 18,2097 
1 1 1 4741 

I9.JM26.. 11.W41 

| S 94444 9 2SS26 

14 62963 6 94444 

I 2 31461 . 4.623S3 

10 639851 2J14H 




Figure 9. Darcy velocity of fracture flow in FEFLOW using simple pressure gradient of 0.307 KPa for comparison with OpenFOAM. 
Image of the model with slices with high hydraulic conductivity used to represent fractures 



TABLEl. VELOCITY AND PRESSURE OF THE FRACTURE FLOWS SIMULATED IN FEFLOW AND OPENFOAM 



Type of Flow 


Velocity 


Pressure 


Velocity/Pressure 


Hydraulic 


Roughness 


Aperture 


Hydraulic 






Gradient 


Gradient 


Conductivity (K) 


Coefficient 




Aperture 


Hagen Poiseuelle 


1.12 


0.307 


3.64 




NA 


0.1 


0.0075 


Manning Strickler 


1.16 


0.307 


3.77 


NA 


50 


0.1 


NA 


Darcy 


1.15 


0.307 


3.74 


37.5 


NA 


0.1 


NA 


High Conductive Layer 


1.15 


0.307 


3.75 


30 


NA 


0.1 


NA 


















Laminar without 


0.15 


0.414 


0.36 


NA 


NA 


0.1 


NA 


surface roughness 
















Laminar with surface 


0.22 


0.417 


0.53 


NA 


NA 


0.1 


NA 


roughness 
















Turbulent without 


1.09 


0.179 


6.09 


NA 


NA 


0.1 


NA 


surface roughness 
















Turbulent with 


1.19 


0.307 


3.88 


NA 


NA 


0.1 


NA 


roughness aperture 
















5 mm 
















Turbulent with 


1.31 


0.734 


1.78 


NA 


NA 


0.1 


NA 


roughness aperture 
















10mm 

















VI. 3dec Modelling 

Rock fractures and flow of fluid media in them have been studied extensively, but the role of the fractures in 
this matter remains unobserved in most cases. The use of 3DEC in modelling the flow of a fluid under a 
pressure gradient has been studied such that a comparison can be made between the flow of the fluid in case 
of using a discrete element approach as opposed to a CFD approach can be lucrative. The flow of fluids in 
fractures modeled in 3DEC assumes a modified cubic law for flow. This approach can be easily incorporated 
into OpenFoam, nevertheless in case of 3DEC, there is a sharp distinction made between the properties of the 
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rock and the fracture. The deformability of the fracture is larger than that compared to the rock; the shear 
resistance of the fracture due to the roughness can be efficiently incorporated into the model. 
A model can be constructed so that the simulation of flow of fluid in fractures incorporates the deformability 
of the fractures as well as the shear resistance due to the roughness property. Although the simulations in this 
paper assign a roughness coefficient to the flow of fluid, this is partly conjectured; actual conditions will have 
to be calibrated using lab and field studies. The incorporation of an accurate simulation of flow will have to 
incorporate rock properties and most importantly the fracture properties though a yielding constitutive model. 
A joint constitutive model can be setup using 3DEC which can be used to represent a more realistic scenario 
of flow of fluids in fractures and joints. 

VII. Conclusions 

Simulation of fractures using simple methods is possible and such as using a highly conductive layer to act as 
a fracture in a flow simulation is not only possible but efficient. However the properties of the slice or layer 
like the hydraulic conductivity, aperture etc. has to be pre -determined such that they may have as a fracture 
network. From the fracture simulations in FEFLOW using a geothermal injecting well, we see that the 
maximum hydraulic head is produced in the case of Darcy flow law used in discrete features. Following this 
is Manning-Strickler, Hagen-Poiseuille and the layer with high conductivity in order of reducing head. The 
analysis of the Peclet number in the figures on the left shows us that the convective heat transport is limited 
to the wells and the area surrounding the wells only. The Peclet numbers are high in the case of Darcy and 
Manning-Strickler flow. The comparison of Darcy flow and Manning Strickler flow in discrete features 
shows us that the Darcy flow is more efficient in terms of computational efficiency and Peclet number. 
The best method to adopt this would be the simulation of fracture networks using a CFD method using a 
simple case and then using this result to preset the properties of the slice or layer in question. 
We have to take cognizance of the fact that the simulation using discrete features to represent a fracture is 
possible as well in FEM codes like FEFLOW. However these processes tend to demand a bit more 
computational power, but nonetheless efficient. The flow itself of the fluid in the fracture is also very 
important and the discrete features help us choose the flow equations which increase the control of the user in 
these simulations. There is a large scope for research in terms of CFD simulation of flow in fractures. The 
determination of the roughness coefficient and the surfaceroughness aperture can be actively studied and 
modeled using FEM as well as CFD to make the simulation of flow in fractures more efficient. 
An alternative presents itself in the case of DEM, the simulation of fractures is efficient with many 
opportunities for modification. The incorporation of such flow into ground water modelling software or 
reservoir simulation software will be beneficial, and this represents scope for future study due to the nature of 
versatility in the application of flow in fractures in 3DEC. Moreover, the flow law used in these fractures 
parallels that used in FEFLOW in some cases. This leads to reason where it would be beneficial to model the 
flow of fractures using DEM. 

References 

[1] W. Riihaak, V. Rath, and C. Clauser, "Detecting thermal anomalies within the Molasse Basin, Southern Germany," 

Hydrogeology Journal, 18, 8, pp. 1897-1915, 2010. 
[2] M. Ilyasov, I. Osterman, and A. Punzi, "Modeling deep geothermal reservoirs: recent advances and future 

problems," Freeden,W., Sonar, T.,Nashed, M.Z. (eds.) Handbook of Geomathematics, p. 679-711, 2010. 
[3] H. Diersch, "Interactive, graphics-based finite-element simulation system - FEFLOW - for groundwater flow and 

contaminant transport processes," Report,WASYGesellschaftfiirwasserwirtschaftlichePlanung und 

SystemforschungmbH, 1994. 

[4] O. Kolditz, and H. Dierch," Quasi-steady-state strategy for numerical simulation of geothermal circulation in hot dry 
rock fractures," International journal of non-linear mechanics, 28 (4), pp. 467-481, 2013. 

[5] R. Nelson, Geologic analysis of naturally fractured reservoirs, Gulf Professional Publishing, 2001. 

[6] T. D. van Golf-Racht, Fundamentals of fractured reservoir engineering. Access Online via Elsevier, 1982. 

[7] J. Bear, C. F. Tsang, and G. De Marsily, Flow and contaminant transport in fractured rock. Access Online via 
Elsevier, 1993. 

[8] O. Kolditz, "Modelling flow and heat transfer in fractured rocks: Conceptual model of a 3-D deterministic fracture 

network," Geothermics, 24(3), 451-470, 1995. 
[9] J. Noorishad, M. S. Ayatollahi, and P. A. Witherspoon, "A finite-element method for coupled stress and fluid flow 

analysis in fractured rock masses," International Journal of Rock Mechanics and Mining Sciences &Geomechanics 

Abstracts, Vol. 19, No. 4, pp. 185-193, August 1982. 



124 



[10] W. Riihaak, P. Schatzl, A. Renz, and H. J. G. Diersch, "Numerical modeling of geothermal processes: Issues and 
examples," 10th International Mine Water Association Congress, Mine Water and the Environment, Karlovy Vary, 
2008. 

[11] C. Clauser, and H. Villinger, "Analysis of conductive and convective heat transfer in a sedimentary basin, 
demonstrated for the Rheingraben," Geophys Journal International, 100 pp. 393-414, 1990. 

[12] S. Sarkar, M. N. Toksoz, and & D. R. Burns, "Fluid flow modeling in fractures," Massachusetts Institute of 
Technology, Earth Resources Laboratory, 2004. 

[13] J. R. Shwechuk, "Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator," Applied 
Computational Geometry: Towards Geometric Engineering (MingC. Lin and Dinesh Manocha, editors), Volume 
1 148 of Lecture Notes in Computer Science, pages 203-222, Springer- Verlag, Berlin, May 1996. 

[14] M. Anderson, "Heat as a ground water tracer," Ground Water, 43 (6), pp. 951-968, 2005. 

[15] H. C. Chen, and V. C. Patel, "Near-wall turbulence models for complex flows including separation," AIAA 
journal, 26(6), 641-648, 1988. 

[16] OpenFOAM® - The Open Source Computational Fluid Dynamics (CFD) Toolbox. 2013. OpenFOAM® - The 
Open Source Computational Fluid Dynamics (CFD) Toolbox. [ONLINE] Available at: http://www.openfoam.com/. 
[Accessed 10 June 2013]. 



125 



